Quantum master equations for non-linear optical response of molecular systems 
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Generalized master equations valid for the third order response of an optically driven multi-level 
electronic system are derived within Zwanzig projection formalism. Each of three time intervals of 
the response function is found to require specific master equation and projection operator. Exact 
cumulant response functions for the harmonic profiles of the potential energy surfaces leading to 
Gaussian spectral diffusion are reproduced. The proposed method accounts for the nonequilibrium 
state of bath at the border of intervals, and can be used to improve calculations of ultrafast non-linear 
spectra of energy transferring systems. 



I. INTRODUCTION 

Non-linear response theory forms the theoretical ba- 
sis for description of many important physical phenom- 
ena and experimental techniques, in particular the multi- 
dimensional techniques of non-linear optical spectroscopy 
In recent years, two-dimensional (2D) coherent spec- 
troscopy , developed first in NMR S , has been brought 
into the near infra-red and optical jj, HJ domains. The 
2D Fourier transformed spectrum completely character- 
izes the third order non-linear response of a molecular 
ensemble in amplitude and phase [2] providing thus, the 
maximal information accessible in a three pulse optical 
experiment. Since its first experimental realization, it has 
yielded new insights into photo- induced dynamics of elec- 
tronic excited states of small molecules [f| , polymers , 
large photosynthetic aggregates H and even solid state 
systems [9]. For instance, long living electronic coherence 
has been detected in energy transfer dynamics of photo- 
synthetic proteins [To|, a find that stimulated a renewed 
interest in the properties of energy transfer in biological 
systems and its arguably quantum nature. 

Most 2D experiments are well described by the third 
order perturbation semiclassical light-matter interaction 
response function theory. Simulations must include de- 
phasing of electronic coherence and electronic energy 
transfer between levels which are both induced by inter- 
action of the electronic degrees of freedom (DOF) with 
environment (bath). Response functions of model few 
level systems with pure dephasing (no energy transfer 
between the levels) can be calculated using second cumu- 
lant in Magnus expansion [lj, and expressed in terms of 
the energy gap correlation function (EGCF), C(t). For 
Gaussian stochastic bath, this result is exact, and the 
EGCF of the electronic transitions thus determines fully 
both the linear and the non-linear (transient) absorption 
spectra. For energy transferring systems, such as Frenkel 
excitons in photosynthetic complexes the exact re- 
sponse functions cannot be constructed by second cumu- 
lant. Photosynthetic complexes are relatively large, and 
the proper methods to simulate finite timescale stochas- 
tic fluctuations at finite temperatures [12| carry a sub- 
stantial numerical cost. Stochastic simulations can thus 
be readily implemented only for small systems fl3l [li) . 



An explicite inclusion of some important environmental 
modes into the Hamiltonian also increases the size of the 
system beyond feasibility. Practical calculations thus re- 
quire some type of reduced dynamics where only elec- 
tronic DOF are treated explicitly. A convenient form 
of the reduced dynamics, master equations, is realized 
by deriving kinetic equations for expectation values of 
relevant quantities averaged over the DOF of the bath 
[l5l floT | . If these quantities (such as transition dipole mo- 
ment) do not depend on the bath DOF, all relevant infor- 
mation is carried by the reduced density matrix (RDM). 

It was demonstrated that RDM master equations de- 
rived by projection operator technique reproduce the lin- 
ear response exactly fl7l |. For higher order response, 
however, the same approach neglects bath correlation 
between different time intervals of the molecular photo- 
induced evolution separated by interaction with laser 
pulses (see e.g. Section 3 in Ref. [Hj). In other words, 
standard RDM master equations do not properly account 
for the non-equilibrium state of the phonon bath present 
after the second and the third laser pulse. This neglect 
can lead to a complete loss of the experimentally observed 
dynamics in simulated 2D spectra, such as in the case of 
the vibrational line shape modulation of a single elec- 
tronic transition [lijj . This modulation interplays with 
simultaneous effect of electronic coherence. Accounting 
reliably for bath correlations among intervals is therefore 
of utmost importance for interpretation of the role of 
quantum coherence in natural light-harvesting. The dif- 
ference between the 2D lineshape calculated with exact 
formula, and by RDM, i. e. neglecting the correlations, 
is demonstrated in Fi g. [T] We employed the definitions 
2D spectra from Ref. [19| and the response theory from 
Ref. Q. 

In this Letter, we show that the exact third order re- 
sponse function of a multilevel molecular system with 
pure dephasing can be calculated by certain generaliza- 
tion of RDM master equations. We introduce special 
parametric projection operators that enable us to derive 
distinct equation of motion for each time interval of the 
response function. These parametric projectors offer a 
systematic approach to improvement of current master 
equation method by incorporating bath correlation ef- 
fects into the calculation of non-linear response in energy 
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FIG. 1: Two dimensional spectra of a two-level electronic 
system interacting with bath of harmonic oscillators at pop- 
ulation time T — 0, left panel: calculated by exact cu- 
mulant Eq. (JSJ, right: by applying RDM for optical co- 
herences (Eq. (|25p ). Frequency- frequency correlation spec- 
trum at T — defined in Ref. [l9l | is shown for line- 
broadening function of overdamped Brownian oscillator [l[ 
g(t) = (AfesTVcorrfr -1 - iAr corr )[e~ t//Tcorr - 1 + t/r COII ] where 
A = 100 cm -1 , T corr = 100 fs, and the temperature is T = 300 
K. 



transfer systems. 



II. MODEL HAMILTONIAN 

We model a molecular aggregate (e.g. photosynthetic 
antenna) as an assembly from N electronic two-level 
molecules coupled by resonance interaction. Third order 
optical experiments on such systems address only the col- 
lective electronic ground state \g) and the states \i) with 
one or two excited molecules of the aggregate. Thus the 
TV-molecule aggregate has M = N+(N-i)x N/2 rele- 
vant excited states. The free evolution of the system be- 
tween the times of interaction with external light pulses 
is described by Frenkel exciton Hamiltonian 

M 

H = Z + V g + e g \g)(g\ + + A^i)^ + H ies . (1) 

t=i 

Here, Z represents kinetic energy of the bath DOF, e g 
and Si are the electronic energies of the ground- and the 
excited states, respectively, A$i = $^ — $ g — ($; — $ g ) is 
the energy gap operator, and $ g and $i are potential en- 
ergy surfaces (PES) in the electronic ground- and excited 
states, respectively. Harmonic PES represents Gaussian 
fluctuations of excitons. The energy gap operator is de- 
fined so that it vanishes (A$j) = trB{A$iW eq } = in 

the ground state equilibrium W eq = e k s T B /zq. Here, 
the ground state bath Hamiltonian reads Hb = Z + <t g . 
The operators A$i describe a macroscopic number of 
DOF so that one cannot easily treat them explicitely in 



numerical calculations. Resonance coupling 

M 

#rcs= Yl J ^G\ ( 2 ) 

between the excited states causes electronic energy trans- 
fer between molecules. Diagonalization of average Hamil- 
tonian (H) is required to identify the positions of peaks 
in absorption spectrum, which correspond to transition 
frequencies w jg = (ej — e g )/h between electronic eigen- 
states \i) with eigenenergies e^, i — 1,...,M and the 
ground state. In the very same spirit, 2D spectra are con- 
veniently interpreted in terms of energy transfer between 
these energetic eigenstates. When resonance coupling be- 
tween molecules in the aggregate vanishes Jy — > the 
eigenstates \i) coincide with excited states \i), and so do 
the eigenenergies and energies £j. 

III. LINEAR RESPONSE 

For harmonic PES and — > the calculation of re- 
sponse functions of all orders can be carried out exactly 
by second order cumulant. The linear response function 

I(t)=Yl\ d i\ 2 Pi9it) + C.C. (3) 

i 

is generated by time evolution of coherence elements of 
the RDM p ig (t) = e -^i S t-gn(t) ^ which can be differenti- 
ated into the convolutionless master equation fl7l |20| 

d 

—p lg = [-iui ig - g u (t)] p ig {t) . (4) 

In Eq. §S§ , di are the transition dipolc moments between 
eigenstate \i) and ground state \g). the line broadening 
function 

9ij{t) = ¥ I dt> I dt " c ^ ^ 

is related to the EGCF, 

C l3 {t) = tr B {e^' h)flB ^ ie - {lt/K)flB AQjWeq}, (6) 

in the ground state equilibrium, and gij = (dgij/dt). 
Further on in this work, we assume that the electronic 
energy gap operators A$i are nonzero. 

IV. MASTER EQUATIONS FOR NON-LINEAR 
RESPONSE 

In non-linear experiment, molecular system is probed 
by three short laser pulses separated by time intervals r, 
T, and its response is monitored after time t. Non-linear 
response function is conveniently dissected into Liouville 
space pathways jlj, each representing one combination 
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FIG. 2: Double sided Feynman diagram of the Liouville path- 
way R 3 2 for a system with three levels \g), \i) and \j). Inter- 
actions with external field are denoted by straight arrows, the 
wavy arrow corresponds to emitted field. The projectors "Po, 
V T and Pyi T are shown explicitly with point of their appli- 



cation in the response function. 



of indices of RDM during the periods r, T, and t. For 
our discussion we select a pathway usually denoted as i?2 
(see Fig. [2] for the corresponding double sided Feynman 
diagram), i.e. the pathway which is responsible for pho- 
ton echo signal, and which contains excited state evolu- 
tion during the second interval (waiting time) T [l| . The 
R2 response function involving excited states \i) and 



|j) reads 



Ri ji \t,T,r) = \d i \ 2 \d J f 



x tr B {W iff (i)W J - i (T)W 9i (r)W eq }. 



(7) 



Here, Uji(t) = Ujijiit) are abbreviated matrix elements 
of the evolution superoperator U(t) which is the Green's 
function solution of the Liouville- von Neumann equation 
with Hamiltonian, Eq. ([1]). Between each two consecu- 
tive Green's functions in the trace of Eq. (j7), the action 
of the external light causes transition between electronic 
levels \g), \i) and For a harmonic PES the second 
order cumulant is exact, and it yields [2l[ 



4 Ji) (t,T,r) 



\di\ a \d s 



-i(u>jgt+{u>jg—U)i g )T-U} ig T) 



x e -9» (r+T)-g n (T+f ) + <;*,■ (t+T+r ) 
Xe -Sy(*)-S«( T )+S , «( T ). (g) 

By differentiating Eq. ^ with respect to t we find the 
response functions to be solutions of evolution (master) 
equations 
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+ IC { i l \t,T,r)}p jg (t;T,T) (9) 



with the relaxation tensor 



K,{ ji 



(t,T,r) = -gU t ) + gUt + T + T)- gjj (t+T). (10) 



and initial condition pj g (0; T,t) = r!£*\q,T,t). Master 
equation can be similarly found in the second interval by 



differentiation of R^ l) (0, T, r) 
d 



pji(T;r) = [-i(ujg - w lg ) + /C^ l) (T, r)] fti (T; r), 



dT 

with relaxation tensor 



(11) 



K,r ) (T,T)=g ij (T)-g* i (T + T) 



(12) 



and initial condition Pj g (0; r) = R 2 (0, 0, r). In the first 
interval we find an analogical master equation for p g i (r) 
which is (up to a complex conjugation) identical to Eq. 
d3|) . We have thus demonstrated that evolution equations 
can be constructed with the structure of the convolution- 
less master equation, and that they yield exact third or- 
der response function. Their forms, however, have now 
to be specific to each interval and to each Liouville space 
pathway. 



V. CONSTRUCTING MASTER EQUATIONS 
BY PARAMETRIC PROJECTORS 

In Ref. [l7j it was pointed out that Eq. (0} agrees 
with convolutionless master equation (CLME) of Hashit- 
sume et al. j22|, when it is evaluated to second or- 
der in A$. Similar microscopic derivation of Eqs. @ 
to (|12|) by using projection techniques would open a 
way to constructing improved master equations for the 
electronic energy transfer in the cases where no exact 
cumulant solution is known. The Zwanzig projection 
technique requires to define projection superoperator by 
its action on a arbitrary operator A. In general form 
V A = ^ ab W a bt r f'B{{a\A\b)}. where the normalized bath 
matrix (irs{W a 6} = 1) can be chosen to achieve specific 
goals. Master equation for the linear response, Eq. (|4]), is 
homogenous and it corresponds to the choice W a b = W eq , 
when QW(0) = (Q = 1 - V), and thus the inhomo- 
geneous part of the CLME vanishes [l7| for uncorrelated 
initial state W(0) = p(0)W eq . In the second and third 
interval the choice W a & = W eq will not yield Eqs. @ 
and (fTTj) , since the inhomogeneous part of the CLME 
does not vanish. In fact, the inhomogeneous terms can- 
not be eliminated in all three intervals simultaneously by 
any single projector. This is a microscopic counterpart of 
the fact that non-linear optical response function cannot 
be exactly evaluated by simulating the reduced dynam- 
ics of the system by a single master equation. Indeed, 
the homogeneous part of the CLME corresponds to only 
one of the four (nonzero) contributions of the third order 
response function 



R^\t,T,r) = tr B {U jg (t) 



■x(V + Q)U 3l {T)(V + Q)U gi ( T )VW eq }. (13) 
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The master equations (|3} to (fT2")l can be derived when we 
allow different projection operator to be defined for each 
time interval of the response function. Let the projection 
operators Vq, V t and V^]_ T in the respective intervals be 
chosen so that they eliminate the inhomogeneous terms in 
the CLME master equation, i.e. so that (1 — Vq)W(0) = 
0; (1-V t )W(t) = 0, etc (see Fig. EJ. The corresponding 
projectors read as follows. The first coherence projection 
superoperator 

V W{t) = W eq tr B {W(t)} = W eqP (t), (14) 

is the widely used Argyres-Kelly projector (23| . The sec- 
ond interval can be treated with the population projection 
superoperator 

V T W{t) =J2pki(t)U gl (r)W eq ^\k){ll (15) 

kl 

which yields a single master equation for all RDM ele- 
ments in the second interval. Projector, Eq. (|15p . to be 
applied to a single RDM clement as in Eq. reads 

Vi ji) Wji{t) = p ]l {t)U gl {T)W eq & ) . In the third interval 
of the response a set of second coherence projection super- 
operators reflecting different bras ((i|) during the second 
interval have to be defined as 

V^ +T W(t) =Y,P^)Uki{T)U gi {r) 

k 

x W eq p^ T \k)(g\. (16) 
The factors (3 { T k) and 

are present to ensure that the 
superoperators have the projection property (V T ) 2 = V T 
and {V% T f — V^\ T . From these conditions it follows 
e.g. p<% = 0® = l/tr B {U gl {T)W eq }. The factor 0® 

evaluates in the second cumulant to /3r ' > = e 9ii ^- T \ 

An alternative to this formalism would be to account 
directly for the Q terms in Eq. (TT51). a nd to use the 
time convolutionless operator identity [22j for the deriva- 
tion of the non-linear response as in Ref . [26] . The two 
approaches lead to the same exact results for a system 
of non-interacting molecules. For a coupled system, the 
quality of the equations of motion derived by the two ap- 
proaches has to be compared in a separate study. How- 
ever, the central idea of the projector approach, i.e. its 
choice in a form which satisfied (1 — V)W = 0, is gen- 
eral, and it does not limit the present method to the 
three projectors suggested above. Depending on physi- 
cal situation, we expect that specialized projectors can 
be introduced for interacting systems, possibly based on 
the comparison with the results of Ref. [26j . 

VI. RECOVERING THE CUMULANT 
EXPRESSIONS 

Let us now demonstrate how these projection operators 
can be used for calculation of the non-linear response 



function, Eq. (J7J), in terms of the reduced dynamics of 
the system. We can write 

R { 2 ^(t,T,T) = \dj\ 2 \di\ 2 

x Ufg (t;T,r )Uji (T; r )U gi (r)p Q , (17) 

where the three reduced evolution superoperators de- 
noted by U are solutions of three different master equa- 
tions following from the projectors, Eqs. (fT4")l to dTHJ) • To 
simplify the matters, we will concentrate on the response 
function where i — j, i.e. on . Since we deal with 
a multi-level electronic system with adiabatically sepa- 
rated states (Jij = 0), we can assume that Uu(T; r) = 1 
and thus 

R^ } (t,T,r) = \di\ 4 a$(t;T,T)hi gi (T)p . (18) 

Note that all T-dependence in Eq. (JT5J is due to the 
second coherence parametric projector. 

Now we construct the master equations for the first and 
second coherence intervals (i.e. in intervals characterized 
by time r and t, respectively), and demonstrate that they 
lead to correct response function. First, we split the to- 
tal Hamiltonian into a sum of the pure bath (denoted 
by B), pure electronic (5) and system-bath interaction 
(S-B) parts. Bath Hamiltonian has already been defined 
above, the electronic Hamiltonian reads H$ = £ g \g){g\ + 
Si e i|*)(*li an d the system-bath coupling Hamiltonian is 
Hs-b = A$|z)(i|. We define the Liouville superop- 
erator Cs-b using the corresponding Hamilton operator 
as £s-b = h 1 [Hs-Bi ■ ■ ■ ]- an d switch to the interac- 
tion picture induced by Hq = Hb + Hs- In the second 
order in £s-b, and assuming that Q x W(0) — for x — 0, 
and i = r + rwe obtain the following quantum master 
equation 

| V x W [ x I] {t) = -tP x Cs-B(t)V x W^(t) 

— J dt'V x Cs-B(t)Q x Cs-B (t')V x W^(t). (19) 
o 

The time-dependence of the Liouvillian and the upper in- 
dex (/) of the density operator denote interaction picture 
with respect to the Hamiltonian Hq = Hb + Hs ■ 

In the first coherence interval we use the projection 
superoperator V§. The construction of the oper- 
ator ensures that trB{A&iW eq } — 0, and so the first 
term in Eq. (|19p contributes with zero. The second 
term is only non-zero when no H$-b occurs on the right 
hand side of the density operator, and thus only the 
term VQH s ^B{t)Hs-B{t')V W {I) (t) contributes. The 
evaluation leads to Eq. ([J} and consequently yields 
^gi( T ) — e" 9ii ^ +l ^ iaT . A straightforward application 
of the second coherence projection superoperator V^ +r 
leads to 

§-p$>(f) = -I T+ r{t)p^ g \t) M T+T (t)p$(t), (20) 
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where 

I T+T (t) = itr B {A^{t)Uu(T)U gl (T)W eq }e^<- T \ (21) 
and 

M T+T (t) = J dt'tr B {A$>(t)A$>(t')Uu{T) 
o 

xU gl (T)W eq }e^ - J dt'I T+T {t)I T+T {t'). (22) 
o 

The coefficients I and M, Eqs. ([HJ and ([2"2)l. can be 
recast into the second cumulant by a technique described 
in Ref. [2"3| . A tedious but straightforward calculation 
gives 

/T+r (t) = <?« (* + T) - g u (t) -g* i {t + T + r)+ gl (t) (23) 
and 

M T+T {t) = gu{t). (24) 

Combining the results of the second order cumulant ex- 
pansion of Eqs. (|2ip and (f2"2")l and switching back from 
the interaction picture we arrive at the master equation 
Eq. (fTUjl . The solution of the master equation can be 
easily obtained in form of the time evolution superoper- 
ator Ui g (t] T, r). Finally, inserting U g i(r) and Ui g (t; T, r) 
into Eq. (fT8f we obtain Eq. (JSj) which is the desired re- 
sult. Calculation can be repeated for any Liouville space 
pathway reproducing the results of direct cumulant ex- 
pansion. 

Note that standard choice of ground state equilibrium 
for projector would yield response function 

Rf\t,T,r) = |d i | 2 |d J | a e- <(Wi « t - Wi « T) 

xe -9U-r)-9ii(t) (25) 

Comparison between 2D lineshapes simulated using Eq. 
© and Eq. (g5j) is made in Figure [TJ 

The projector technique developed above can be ex- 
tended to finite resonance coupling by straightforward 
abandoning of the — > limit. We will show elsewhere 
(25| that provided the optical coherences in the first inter- 
val can be assumed independent of each other (secular ap- 
proximation), the projectors, Eqs. (|14j) and (|15|) can be 



used to derive corresponding master equations. The com- 
plete response function involves a single master equation 
for the first coherence interval and two master equations 
for the waiting interval. In the second coherence interval 
of the response function there are iV master equations 
in each Liouville pathways, where N is the number of 
single exciton states. When energy transfer between sys- 
tem's eigenstates is induced by non-zero couplings, 
the application of the parametric projectors represents 
an approximation whose quality depends on the previ- 
ous energy transfer history of the system. Nevertheless, 
unlike the usual master equations based on standard pro- 
jectors, the parametric projectors lead to correct limit in 
Jij — > 0. The prescription (1 — V)W(t) — for construc- 
tion of the parametric projection operators is general, 
and can be used to construct a proper projector in the 
cases where V t +t would fail. 

VII. CONCLUSIONS 



In this Letter, we have outlined a general method for 
calculation of non-linear response functions by quantum 
master equations. We have tailored parametric projec- 
tion superoperators to derive master equations for each 
time interval of the response function, and we demon- 
strated the validity of the method by reproducing a well- 
know result exact for a multi-level system. The results 
presented here can be extended to a master equation 
method for calculation of a general multi-point correla- 
tion functions. This opens a new research avenue with 
consequences in the theory of open systems, non-linear 
spectroscopy and potentially other branches of physics. 
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